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Abstract 

This paper aims to compare rational Chebyshev (RC) and Hermite functions (HF) collo- 
cation approach to solve the Volterra's model for population growth of a species within a 
closed system. This model is a nonlinear integro-differential equation where the integral 
term represents the effect of toxin. This approach is based on orthogonal functions which 
will be defined. The collocation method reduces the solution of this problem to the solu- 
tion of a system of algebraic equations. We also compare these methods with some other 
numerical results and show that the present approach is applicable for solving nonlinear 
integro-differential equations. 
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1. Introduction 

Many science and engineering problems arise in unbounded domains. Different spectral 
methods have been proposed for solving problems in unbounded domains. The most common 
method is through the use of polynomials that are orthogonal over unbounded domains, such 
as the Hermite spectral method and the Laguerre spectral method [H- 



Guo [1CH12| proposed a method that proceeds by mapping the original problem in an 
unbounded domain to a problem in a bounded domain, and then using suitable Jacobi 
polynomials to approximate the resulting problems. 

Another approach is replacing infinite domain with [— L, L] and semi-infinite interval 



with [0, L] by choosing L, sufficiently large. This method is named domain truncation 13 
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Another effective direct approach for solving such problems is based on rational approx- 



imations. Christov [14| and Boyd [15|, |16| developed some spectral methods on unbounded 
intervals by using mutually orthogonal systems of rational functions. Boyd [l6[ defined a 
new spectral basis, named rational Chebyshev functions on the semi-infinite interval, by 
mapping to the Chebyshev polynomials. Guo et al. [l7| introduced a new set of rational 
Legendre functions which are mutually orthogonal in L 2 (0, +00). They applied a spectral 
scheme using the rational Legendre functions for solving the Korteweg-de Vries equation on 



the half line. Boyd et al. [18| applied pseudospectral methods on a semi-infinite interval 
and compared rational Chebyshev, Laguerre and mapped Fourier sine. 

Parand et al. 19H25| applied spectral method to solve nonlinear ordinary differential 
equations on semi-infinite intervals. Their approach was based on rational Tau and colloca- 
tion method. 

Among these, an approach consists of using the collocation method or the pseudospectral 
method based on the nodes of Gauss formulas related to unbounded intervals @. 

Collocation method has become increasingly popular for solving differential equations 
also they are very useful in providing highly accurate solutions to differential equations. 

We aim to compare rational Chebyshev collocation (RCC) approach and Hermite func- 
tions collocation (HF) collocation approach to solve a population growth of a species within 
a closed system. 

This paper is arranged as follows: in subsection II. H the Volterra's population model is 
taken into consideration and some of the traditional methods that solved it are discussed. 
In sections [2] and [31 we describe the properties of rational Chebyshev and Hermite functions. 
In section HJ the proposed methods are applied to solve Volterra's population model. This 
equation is first converted to an equivalent nonlinear ordinary differential equation and then 
our methods are applied to solve this new equation, and then a comparison is made with 
existing methods that were reported in the literature. The numerical results and advantages 
of the methods are discussed in the final section. 

1.1. Volterra's Population Model 

The Volterra's model for population growth of a species within a closed system is given 
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as 



— = ap — bp — cp / p(x)dx, p(0) = po, (1) 
at J 

where a > is the birth rate coefficient, b > is the crowding coefficient, and c > is 
the toxicity coefficient. The coefficient c indicates the essential behavior of the population 
evolution before its level falls to zero in the long term. p is the initial population, and 
p = p{t) denotes the population at time t. 

This model is a first-order integro-ordinary differential equation where the term cp J Q t p(x)dx 



2 



represents the effect of toxin accumulation on the species. We apply scale time and popula- 
tion by introducing the nondimensional variables 



tc pb 
t — — and u = — , 
b a 

to obtain the nondimensional problem 
du 

k— = u — u 2 — u I u{x)dx, u(0) = uq, (2) 
dt Jq 

where u(t) is the scaled population of identical individuals at time t, and k = c/(ab) is a 
prescribed nondimensional parameter. The only equilibrium solution of Eq. fl2]) is the trivial 
solution u(t) — and the analytical solution [28| 

u(t) = uq exp(— I [1 — u(t) — I u(x)dx]dr), 
K Jo Jo 

shows that u(t) > for all t if u > 0. 

The solution of Eq. ([!]) has been of considerable concern. Although a closed form solution 



has been achieved in |26l . |27| | , it was formally shown that the closed form solution cannot 



lead to any insight into the behavior of the population evolution 26]. In the literature, 
several numerical solutions for Volterra's population model have been reported. In 26], 
the successive approximations method was suggested for the solution of Eq. (|2J), but was 
not implemented. In this case, the solution u(t) has a smaller amplitude compared to the 
amplitude of u(t) for the case /t C 1. 

In [27j], the singular perturbation method for solving Volterra's population model is 
considered. The author scaled out the parameters of Eq. as much as possible and 
considered four different ways to do this. He considered two cases k = c/(ab) small and 
K — c/ (ab) large. 



It is shown in [27[ that for the case ft < 1, where populations are weakly sensitive to 
toxins, a rapid rise occurs along the logistic curve that will reach a peak and then is followed 
by a slow exponential decay. And, for large k, where populations are strongly sensitive to 
toxins, the solutions are proportional to sech 2 (t). 

In [28( , several numerical algorithms namely Euler method, modified Euler method, clas- 
sical fourth-order Runge-Kutta method and Runge-Kutta-Fehlberg method for the solution 



of Eq. (j2J) are obtained. Moreover, a phase-plane analysis is implemented. In [28|, the nu- 
merical results are correlated to give insight on the problem and its solution without using 
perturbation techniques. However, the performance of the traditional numerical techniques 
is well known in that it provides grid points only, and in addition, it requires large amounts 
of calculations. 

In |29(] Adomian decomposition method and Sinc-Galerkin method were compared for 
the solution of some mathematical population growth models. 
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In 30] , the series solution method and the decomposition method are implemented inde- 
pendently to Eq. (EJ) and to a related nonlinear ordinary differential equation. Furthermore, 
the Pade approximations are used in the analysis to capture the essential behavior of the 
populations u(t) of identical individuals and approximation of u max and the exact value of 
u m ax for different k were compared. 

The authors of 19N21| applied spectral method to solve Volterra's population on a semi- 
infinite interval. This approach is based on a rational Tau method. They obtained the 
operational matrices of derivative and the product of rational Chebyshev and Legendre 
functions and then applied these matrices together with the Tau method to reduce the 
solution of this problem to the solution of a system of algebraic equations. 

In [3l| second derivative multistep methods (denoted SDMM) are used to solve Volterra's 
model. They first converted the model to a nonlinear ordinary differential equation and then 
the new SDMM applied to solve this equation. 



In 32 



the approach is based upon composite spectral functions approximations. The 
properties of composite spectral functions consisting of few terms of orthogonal functions 
are utilized to reduce the solution of the Volterra's model to the solution of a system of 
algebraic equations. 



In 33 



a numerical method based on hybrid function approximations was proposed 

to solve Volterra's population model. These hybrid functions consist of block-pulse and 
Lagrange-interpolating polynomials. 
Momani et al. 
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used a numerical and Analytical algorithm for approxi- 
mate solutions of a fractional population growth model respectively. The first algorithm is 
based on Adomian decomposition method (ADM) with Pade approximants and the second 
algorithm is based on homotopy analysis method (HAM). 

In total, in recent years, numerous works have been focusing on the development of more 
advanced and efficient methods for initial value problems especially for stiff systems. 



2. Rational Chebyshev Functions 



This section is devoted to introducing rational Chebyshev functions (which we denote 
(RC)) and expressing some basic properties of them that will be used to construct the RC 
collocation (RCC) method, rational Chebyshev functions denoted by R n (x) are generated 
from well known Chebyshev polynomials by using the algebraic mapping 4>(x) = (x — L)/(x + 



l) 0, H, 0, M 



R n (x) = T n ((f)(x)), 



(3) 



where L is a constant parameter and T n (y) is the Chebyshev polynomial of degree n. The 



constant parameter L sets the length scale of the mapping. Boyd [l_3j, |37| offered guide- 
lines for optimizing the map parameter L where L > 0. Using properties of Chebyshev 
polynomials and RC, we have 
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Other properties of RC and a complete discussion on approximating functions by RC are 



given m 



22, 36 



2.1. Rational Chebyshev functions approximation 
Let 

at = span{R , R 1: Rn}- 
We define P N : L 2 W (A) 3l N by 

AT 

P N u(x) = y^a fc -Rfc(^) 

fc=0 



(5) 



(6) 



To obtain the order of convergence of rational Chebyshev approximation, we define the space 



H^ A (X) — {v : v is measurable and |H| r ,z,A < oo}, 
where the norm is induced by 



(7) 



\ v Wr.x.A 



E 



(x + l)i 



. d k 

dx k 



2 \ 2 



and A is the Sturm-Liouville operator as follows: 

Av{x) = —w~ l (x)-^- ( w^ 1 (x)— v (x) ) . 

ax V a.x / 



(9) 



is the weight function and w(x) = a/L/(\/x(x + L)). We have the following theorem 
for the convergence: 



Theorem 1. For any v G A (A) and r > 0, 

\\ P NV - v\\ w < cN~ r \\v \\ r ,w,A- 

Proof 1. A complete proof is given by Guo et al.jsdl]. 



(10) 



This theorem shows that the rational Chebyshev approximation has exponential conver- 
gence. 
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3. Properties of Hermite Functions 

In this section, we detail the properties of the Hermite functions (HF) that will be used to 
construct the Hermite functions collocation (HFC) method. First we note that the Hermite 
polynomials are generally not suitable in practice due to their wild asymptotic behavior at 
infinities jjjl ]. 

Hermite polynomials with large n can be written in direct formula as follow: 

„ , , r(n+l) x 2 /2 . r- nix 

~ n"' 2 e''' 2 cos(V2nTTa: - 

Hence, we shall consider the so called Hermite functions. The normalized Hermite functions 
of degree n is defined by 

~ 1 2 

HJx) = - e~ x /2 HJx), n>0,iGl. 

Clearly, {H n } is an orthogonal system in L 2 (M),i.e., 



H n (x)H m (x)dx = a/7t5„ 



where 5 nm is the Kronecker delta function. 

In contrast to the Hermite polynomials, the Hermite functions are well behaved with the 
decay property: 

\H n (x) \ — > 0, as \x\ — > oo, 
and the asymptotic formula with large n is 

H n (x) ~ n~i cos(\/2n + lx — ) 

The three-term recurrence relation of Hermite polynomials implies 



H n+l (x) = x^J —^H n (x) - <f^Hn-i(x), n > 1, 
H (x) = e~ x2 / 2 , Hx(x) = V2xe~ x2 / 2 . 
Using recurrence relation of Hermite polynomials and the above formula leads to 

H' n (x) = \f 7 lnH n _ 1 {x) - xH n (x) 

Ti — / Tl ~\~ 1 ~ 

-H n -i(x) - d -^—H n+1 (x). 



6 



and this implies 



v/ n(n— 1)tt 
2 



H' n {x)H' m {x)dx = < 



m = n — 2, 
0r(n+ 1), 111 = 71, 

m = n + 2, 
0, otherwise. 



2 - 



Let us define 



P N ■= { u : u = e~ x2/2 v,\fv G Pjv}. 



where P/v is the set of all Hermite polynomials of degree at most N. 

We now introduce the Gauss quadrature associated with the Hermite functions approach. 

Let {xj}^ =0 be the Hermite-Gauss nodes and define the weights 



w ; 



71 



Then we have 



(N + l)H%( Xj ) 



N 



0<j<N. 



For a more detailed discussion of these early developments, see the [39|, |40 



3.1. Approximations by Hermite Functions 
Let us define A := {x\ — oo < x < oo} and 

Hn = span{H (x), H%(x), H N (x)} 

The L 2 (A)-orthogonal projection £ N : L 2 (A) — > 1-L N is a mapping in a way that for any 
v e L 2 (A), 



< £nv -v,(f)>= 0, 



G H 



N 



or equivalently, 



iV 



£ N v(x) = y^ViH; 



I (X l 



1=0 



To obtain the convergence rate of Hermite functions we define the space H r A {K) defined by 
H\(K) = {v\v is measurable on A and \\v \\ r A < oo}, 



and equipped with the norm ||u|| r>J A = ||^4 r t>||. For any r > 0, the space H A (A) and its norm 
are defined by space interpolation. By induction, for any non-negative integer r, 

r 

A r v{x) = ^{x 2 + l) {r - k)/2 Pk {x)d k x v{x), 

k=0 

where pk{x) are certain rational functions which are bounded uniformly on A. Thus, 

\\v\\ r ,A <c(j2\\ (X 2 + lt- k ^ 2 Pk (x)d k x V \\\ . 

\fe=0 / 

Theorem 2. For any v e H A (A), r > 1 and < fi < r, 

Proof. A complete proof is given by Guo et al. j4l|. Also same theorem has been proved by 



Shen et al. 38 



4. Solving Volterra's Population Model 

In this section, we study an algorithm for solving Volterra's population model by using the 
collocation method based on rational Chebyshev and Hermite functions. We first convert 
Volterra's population model in Eq. §2§ to an equivalent nonlinear ordinary differential 
equation. Let 

y (t) = j u (x)dx. (12) 
Jo 

This leads to 

y>(t)=u(t),y"(t)=u'(t). (13) 
Inserting Eq. f|T2|) and Eq. f|T3|) into Eq. (T5]) yields the nonlinear differential equation 

Ky"(t)=y>(t)-(y>(t)) 2 -y(t)y>(t), (14) 
with the initial conditions 

2/(0) =0, (15) 
2/(0) =u , 

that were obtained by using Eq. f|T2|) and Eq. f|T3|) respectively. 

We are going to solve this model for u = 0.1 and various n = 0.02, 0.04, 0.1, 0.2 and 0.5. 



8 



4-1. Solving Volterra's Population Model by Rational Chebyshev Functions 

In the first step of our analysis, we apply Pn operator on the function y(t) as follows: 

N 

P N y{t) = J2 a k R k(t) (16) 

fc=0 

Then, we construct the residual function by substituting y(t) by Pjv?/(t) in the volterra's 
population in Eq. 



Res(t) =K^P N y(t) - j t P N y(t) + (j t p Ny{t) 

+ (P N y(t)) (iPNv(t)j , (17) 

The equations for obtaining the coefficients a^s come from equalizing Res(t) to zero at 
rational Chebyshev- Gauss-Radau points plus two boundary conditions: 

(Res( Xj ) = 0, j = l,2,...,iV-l, 

J P N y(0) = 0, (18) 

{£PNy(t)\ t=o = u . 

Solving the set of equations we have the approximating function P^y(x). 



4-2. Solving Volterra's Population Model by Hermite functions 

As mentioned before, Volterra's Population Model is defined on the interval (0, +oo); but 
we know the properties of Hermite functions are derived in the infinite domain (— oo, +oo). 
Also we know approximations can be constructed for infinite, semi-infinite and finite inter- 
vals. One of the approaches to construct approximations on the interval (0, +oo) which is 
used in this paper, is the use of mapping, that is a change of variable of the form 

u = <f)(z) = ln(sinh(foz)). (19) 

where k is a constant. 

The basis functions on (0, +oo) are taken to be the transformed Hermite functions, 

H n (x) = H n (x)o ( j ) (x)=H n ( ( f ) (x)). 
where H n (x) o <fi(x) is defined by H n (x)((p(x)). The inverse map of u = 4>(z) is 



■ln(e w + Ve 2w + 1). (20) 



k 

Thus we may define the inverse images of the spaced nodes {xj} x JZ[_^ as 

T = {^{t) : -oo < t < +00} = (0, +00) 
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and 



x j = (f>- 1 (xj) = e kx ^ j =0,1, 2,... 



Let w(x) denotes a non-negative, integrable, real-valued function over the interval T. We 
define 

L W (T) = {v : T — > M. | v is measurable and || v \\ w < 00} 

where 



00 



I v(x) I w(x)dx 
is the norm induced by the inner product of the space L^(T), 

POO 

<u,v> w = / u(x)v(x)w(x)dx. (21) 
Jo 

Thus {if n (x)} ne N denotes a system which is mutually orthogonal under (J2"lj) . i.e., 

< H n (x),H m (x) > W ( X )= \/n5 nm , 

where = coth(x) and 5 nm is the Kronecker delta function. This system is complete in 
L^(r). For any function / G ^(T) the following expansion holds 



+N 
k=-N 

with 

< /(a;), 

\\ H k {x)\\l (x) 

Now we can define an orthogonal projection based on transformed Hermite functions as 

below: 

Let 

H N = span{H (x), H^x),..., H n (x)} 
The L 2 (r)-orthogonal projection £jv : L 2 (T) — > H n is a mapping in a way that for any 

v e ^ 2 (r), 

< Ijv2/ -y,<p>= 0, V0 G Hat 
10 



or equivalently, 

N 

i N y(x) = J2^H l (x). (22) 

8=0 

Then to apply the HFC method to approximate the Volterra's Population Eq. (|14p with 
initial conditions ( TT5]) . we use the operator introduced in Eq.f l2"2"j) basically. 

We note that the Hermite functions are not differentiable at the point x = 0, therefore 
to approximate the solution of Eq. (fl4|) with the initial conditions Eq. (Tl5|) we construct a 
polynomial pit) that satisfy y'(0) = uq in Eq. (fl5l) and also multiply operator Eq. (|22i) by 
t to make it differentiable at point t = in Eq. (fl5|) . This polynomial is given by 

pit) = \t 2 + 0.1t, (23) 

where A is constant to be determined. 



Therefore, the approximate solution of yit), in Eq. ffT4"l) with initial conditions Eq. ( TT5|) 
is represented by 

? N y(t)=p(t) + ti N y(t), (24) 

Now we construct the residual function by substituting yit) by £,Ny(t) in the Volterra's 
population Eq. f[T4"l) : 

ites,(*) = K ^ N y(t/l) - ^<Ny{t/l) + (J^NV(t/l) 

+ (&v(f/0) (^&w(*/0) . (25) 

where Z is a constant that is called domain scaling. 

It has already been mentioned in 42| that when using a spectral approach on the whole 
real line R one can possibly increase the accuracy of the computation by a suitable scaling of 
the underlying time variable t. For example, if y denotes a solution of the ordinary differential 
equation, then the rescaled function is yit) = y(t/l), where I is constant. Domain scaling is 
used in several of the applications presented in the next section. For more detail we refer 
the reader to (43| . 

The equations for obtaining the coefficients a^s come from equalizing Resit) to zero at 
transformed Hermite-Gauss points: 

Resixj)=0, j = 0, ...,JV + 1, (26) 

where the XjS are N + 2 transformed Hermite-Gauss nodes by Eq. (f20|) . This generates a set 
of N + 2 nonlinear equations that can be solved by Newton method for unknown coefficients 
a,-s and A. 
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Table [T] Shows a comparison of methods in [3jJ, |32j , and the present methods with the 
exact values 



Umax = 1 + « In — (27) 

1 + K — Uo 



evaluated in 28 



Figures [T] and H] show the results of rational Chebyshev and Hermite functions collocation 
methods for k— 0.02, 0.04, 0.1, 0.2, 0.5. These figures show the rapid rise along the logistic 
curve followed by the slow exponential decay after reaching the maximum point and when 
k increases, the amplitude of u(t) decreases whereas the exponential decay increases. 

Figure [3] illustrates a comparison between the two presented methods for k = 0.02. 

Logarithmic graphs of absolute coefficients \ai\ of rational Chebyshev and Hermite func- 
tions in the approximate solutions for k = 0.5 are shown in Figures H] and \5\ respectively. 
These graphs illustrate that both of methods have an appropriate convergence rate. 



5. Conclusions 

The aim of the this study is to develop an efficient and accurate numerical method 
based on orthogonal functions for solving the Volterra model for the population growth of 
a species in a closed system. The methods were used in a direct way on a semi-infinite 
domain without using linearization, perturbation or restrictive assumptions. In this paper, 
we have applied both the rational Chebyshev and the Hermite functions in solving nonlinear 
integro-differential equations and compared the results obtained by the two methods and 



others reported in |31l . |32| . The study showed that rational Chebyshev collocation method 
is simple and easy to use. It also minimizes the computational results. In total an important 
concern of spectral methods is the choice of basis functions; the basis functions have three 
properties: easy to computation, rapid convergence and completeness, which means that 
any solution can be represented. The stability and convergence of rational Chebyshev and 
Hermite functions approximations make this approach very attractive and contributing to 
the good agreement between the approximate and exact values for u max in the numerical 
example. 
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Table 1: A comparison of methods in 31, 32] and the present methods with the exact values for u n 



Present methods 



Other methods 



K 


Exact U max 


N 


HFC 


N 


RCC 


SDMM [31 


CSF with N = 100 [32] 


0.02 


0.92342717 


20 


0.92342704 


14 


0.92342715 


0.92342714 


0.9234262 


0.04 


0.87371998 


25 


0.87371998 


14 


0.87371998 


0.87381998 


0.8737192 


0.1 


0.76974149 


20 


0.76974149 


14 


0.76974149 


0.76974140 


0.7697409 


0.2 


0.65905038 


25 


0.65905038 


11 


0.65905038 


0.65905037 


0.6590497 


0.5 


0.48519030 


30 


0.48519030 


13 


0.48519030 


0.48519029 


0.4851898 
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Figure 1: The results of rational Chebyshev collocation method calculation for k = 0.02,0.04,0.1,0.2,0.5, 
in the order of height 
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Figure 2: The results of Hermite function collocation method calculation for k — 0.02,0.04,0.1,0.2,0.5, in 
the order of height 



16 



1 -1 




HFC <- RCC 



Figure 3: The results of the comparison of rational Chebyshev and Hcrmitc functions collocation method 
for k = 0.02 
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Figure 4: The Logarithmic graph of absolute coefficients |dj| of rational Chebyshev functions in the approx- 
imate solution for k = 0.5 
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Figure 5: The Logarithmic graph of absolute coefficients of Hermite functions in the approximate solution 
for k = 0.5 
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